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Phase-ordering dynamics in nematic liquid crystals has been the subject of much active inves- 
tigation in recent years in theory, experiments and simulations. With a rapid quench from the 
isotropic to nematic phase a large number of topological defects are formed and dominate the sub- 
fi^ . sequent equilibration process. We present here the results of a molecular dynamics simulation of 

Q ' the Gay-Berne model of liquid crystals after such a quench in a system with 65536 molecules. Twist 

C/3 ^ disclination lines as well as type-1 lines and monopoles were observed. Evidence of dynamical scaling 

was found in the behavior of the spatial correlation function and the density of disclination lines. 
However, the behavior of the structure factor provides a more sensitive measure of scaling, and we 
observed a crossover from a defect dominated regime at small values of the wavevector to a thermal 
fluctuation dominated regime at large wavevector. 



C^ 



O 

o 






I 



61.30.Jf, 64.70.Md, 61.30.Cz 



I. INTRODUCTION 
> 

o 

psj ' Topological defects formed during quenches from high-temperature equilibrium phases are of interest in a wide 

T^ [ variety of fields from condensed matter physics to cosmology |l]^]. Uniaxial nematic liquid crystals are excellent 
^O ' materials for studying topological defects because of the variety of defects they possess and because of the ease with 
^D I which they can be studied experimentally. Tables of processes involving defects, such as found in ||5|] are interesting 
^^ r. both from a theoretical and from an experimental point of view. Simulations in which actual molecular configurations 
can be viewed and tracked could greatly elucidate these processes and aid our general understanding of defect dynamics 
and phase ordering. This paper represents a step towards these goals. 

Simulations of defects in nematics have often used, by analogy with 0{n) and other model simulations |^-^, a cell- 
dynamical scheme [prpl in which the order parameter tp at each site is advanced in time according to a time-dependent 
i-rt Ginsburg-Landau equation, within the one elastic constant approximation. Others [0,0| have performed Monte Garlo 

P5 I simulations of a discretized Frank free energy, including allowance for elastic anisotropy and surface anchoring. Still 
O ■ others |p^-p^ have investigated specific types of defects or processes by directly creating the appropriate configurations 
^ ' as initial conditions and then evolving the system. While all of these approaches have yielded fruitful results, it would 
certainly be advantageous to study defects using more realistic off-lattice models with no prior bias towards forming 
J^ ■ any particular defect configurations. In this paper wc present results of a simulation of a quench of the Gay- 
rS I Berne nematic liquid crystal ll^ , a phenomenological fluid model which mimics the behavior of ellipsoidal molecules 
j^ ■ interacting through a combination of attractive and repulsive forces. This model has proven over the past decade 
to capture the essential physical features of real liquid crystals |18{ , and it is an appropriate model for studying the 
formation of topological defects with an off-lattice model. 

This paper is organized as follows. In the next section we review the classification of nematic defects, the dynamical 
scaling hypothesis and the scaling forms of the real-space correlation function and structure factor. In Sec. Ill 
we present the computational details of our simulation, followed in Sec. IV by a description of our defect-finding 
algorithms. Our results and a comparison with theoretical predictions is presented in Sec. V, which is followed in the 
last section by some concluding remarks. 
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II. THEORETICAL BACKGROUND 



A basic understanding of the defects in nematics goes back as far as the early work of Lehmann |19] , but the first 
quantitative classification was given by Oseen [Q. Topological defect solutions are local minima of the Frank free 
energy 
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where fi is the nematic director. A three-dimensional uniaxial nematic has stable point (monopole) and line (discli- 
nation) defects . The former include both radial (charge +1) and hyperbolic (charge -1) geometries pl| which are 
topologically equivalent. Disclination line defects are either of the wedge or twist variety. Either variety is charac- 
terized by a ±180° director rotation about the line (i.e. the defects have charge ±1/2). A twist disclination loop is 
shown in Fig. ^. Note that the loop carries zero monopole charge and the director configuration is uniform at large 
distances from the loop. Wedge disclination loops on the other hand carry a net charge, and at large distances from 
them the director configuration is equivalent to that of a monopole of charge 1. 

Another type of line defect is characterized by director rotations of ±360° (type-1 line)and is unstable to "escaping 
in the third dimension," into a non-singular configuration p2| (see Fig. 0). These escaped structures can still be 
observed experimentally 123-E3], however, and we can also visualize them in our simulation. 
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FIG. 1. (a) Director configuration around a twist disclination line (pointing out of the page). On the left of the defect, the 
director is pointing out of the page, parallel to the disclination line, while on the right, the director lies in the plane of the page. 
Above and below the defect, the directors are depicted in intermediate orientations, (b) Illustration of the difference in director 
orientation between regions interior and exterior to the twist disclination line. The two regions have uniform orientations but 
are rotated with respect to each other by 90° along an axis perpendicular to the loop. 




FIG. 2. (a) Side view of a singular but unstable line with topological charge -1-1 (type-lline). (b) Escaped type-lline with 
directors tilted into a non-singular configuration. 



The dynamical scaling hypothesis |g^ for phase ordering processes asserts that there is a characteristic length L{t) 
{e.g., the domain size or defect separation) such that the system ap pea rs to be time-independent (in a statistical 
sense) when all lengths are rescaled by L{t). For nematics, theory Ig^l predicts that Lit) ~ t^/^, where t is the 
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time since the ordering process began {e.g., the time since a temperature quench which leads to an isotropic-nematic 
transition). The dischnation hue density pdisc (total length of disclination lines per unit volume) should then scale 
as L{t)/{L(t))^ ^ {L{t))^'^ r^ t^^ and the monopole density Pmonop (number of monopoles per unit volume) as 
{L{t))~^ ^ t~^^^ . Note that defects occur at the intersections of domains growing with differing director orientations 
(the Kibble mechanism |g^). Until the domains are large enough, defects are neither well-defined nor well-separated 
(one needs a defect separation larger than the defect core size 11^) so scaling is only assumed to hold at "late" times. 
Experiments are generally consistent with scaling predictions except perhaps in the behavior of monopoles P,Eq,E9[ . 
Simulations also demonstrate scaling [p 10 13] but with calculated exponents somewhat different from theory and 



experiment. There have also been indications in simulations that more than one characteristic length may be present 
|lO| , but this is possibly just a finite-size effect. 

The real-space order parameter correlation function C(r, t) and its Fourier transform ^(k, t) are widely used probes 
|p6| of domain structure and dynamical scaling. For the nematic order parameter 
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one has the definitions 
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According to the scaling hypothesis [^, the data for the orientationally averaged C{r,t) at different times should 
collapse to a single curve when distances at time t are rescaled by L{t). Similarly, S{k,t) should have a single, 
underlying scaling form: that is, 

C{r,t)^f[r/L{t)], (4a) 

S{k,t)^L'^g[kL(t)]. (4b) 

The late-time behavior of S{k) is determined by the type and number of defects present in the system. For nematics, 
S{k) can be written in the form p^] 

_ 367r'' Stt^ 3 fc^T 

'=>\K,t) - Pyno7iop—j^ + Pdisc-j^ + ^J^, (5) 

where the right-hand side includes contributions from the monopoles, disclinations and thermal fluctuations respec- 
tively (the nonsingular type-1 lines do not make a power-law contribution to the structure factor). Here K is the 
elastic constant in the one-constant approximation. For thin nematic films (i.e. two spatial dimensions, but with 
a three-dimensional director n) the monopole and disclination contributions to eqn. (H) are replaced by a single 
contribution proportional to A;~^ arising from disclination points characterized by ±180° director rotations. In the 
three-dimensional case the disclination contribution proportional to k~^ appears for twist and wedge disclination 
loops (as well as any curved disclination loop segment) at wavevectors k ^ R~^ , where R is the radius of the loop 
|30| . For smaller wavevectors, there is no power-law contribution to the structure factor from the twist loops (recall 
that the director configuration is homogeneous at large distances from the loop), and wedge loops contribute a term 
of the same form as the monopoles. The defect contributions to the structure factor above are specific examples of 
Porod's Law pq] which states that 

where p is the defect density, d is the number of spatial dimensions and D is the defect dimensionality {e.g., points have 
D — Q and lines have D — 1). Experiments pl| show good scaling of S{k) with an asymptotic exponent approximately 
equal to 5, with the approach through effective exponents lying between 5 and 6 p^]. Zapotocky and Goldbart |3C| ] 
have shown that such behavior would be consistent with the presence of sufficient numbers of monopoles or wedge 
disclination loops. However, experimentally the population of monopoles seems too low, and wedge disclinations are 
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energetically less preferable [^3|j34[| than twist disclinations for typical values of the nematic elastic constants (though 
the former defects might be generated dynamically). From Eq. (|^) we see that thermal fluctuations will dominate 
the structure factor for sufficiently large wavevectors satisfying {^ksT / K pdisc)k^ > 1 (assuming that monopoles are 
not present in large numbers). This behavior has not been seen in the experimental studies ||35|,^ carried out thus 
far. As discussed in |3Q], the scattering experiments were performed over a time range where the defect density is 
sufficiently large that the crossover to thermal regime is not evident for wavevectors in the visible range. 



III. SIMULATION DETAILS 



We performed a molecular dynamics (MD) simulation using the Gay-Berne model [ [l7[ , an intermolecular potential 
similar to the simple Lennard- Jones potential but extended to model the anisotropic mesogen shape. The complete 
Gay-Berne potential is as follows [ p7[ : 

f/(ui,Uj,f) = 4e(ui,Uj,f) 
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where Ui, Uj give the orientations of the long axes of molecules i and j, respectively, and r is the intermolecular vector 
(r = Ti — Tj). The parameter a (ui, Uj, f) is the intermolecular separation at which the potential vanishes, and thus 
represents the shape of the molecules. Its explicit form is 
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I-X(ui-Uj) 
where ctq — Cs (defined below) and x is 
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Here tTg is the separation between two molecules when they are oriented end-to-end, and Og the separation when 
side-by-side. The well depth e (ui, Uj, f ), representing the anisotropy of the attractive interactions, is written as 



e (ui, Uj, f) = eoe" (ui, Uj) e'^ (ui, Uj, f) 



where 



and 



e(ui,Uj) = {l-x^(ui-Uj)^} 
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with x' defined in terms of Eg and £s^ the end-to-end and side-by-side well depths, respectively, as 
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The overall energy scale is set by the value of eq. Some representative plots of the Gay-Berne potential energy curves 
are shown in Fig. g. 
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FIG. 3. Potential energy curves for the Gay-Berne model using the parameters cited in the text. Curves are shown for four 
sample molecular pair configurations as indicated. 

For the adjustable parameters, we used the values suggested by Berardi et al. [ p8[ : fi — 1,^ — 3,cre/(Js — 3, 
and Es/Se — 5. These values yield a nematic phase over a wider range of temperatures, compared to the original 
parameterization chosen by Gay and Berne ||l^. The temperature was controlled by velocity rescaling [^, and 
the density was fixed at p* — 0.3 with the dimensions of the simulation box in the ratio 2:2:1. Periodic boundary 
conditions were applied. The system was equilibrated at T* = 3.6 in the isotropic phase for 130000 MD time steps 
(with a dimensionless time step of 0.004) and then a quench to T* — 3.2 was implemented (the nematic-isotropic 
transition temperature is approximately 3.5). The system gradually came to equilibrium in the nematic phase with 
the order parameter S (the largest eigenvalue of Qaf^) saturating at a value of 0.69 over the next 100000 steps. We 
used a domain decomposition approach on the Cray T3E at the San Diego Supercomputing Center. Briefly, the 
domain decomposition approach involves dividing up the simulation volume into a number of cells, each controlled by 
a different processor (we used 64 cells). Because the Gay-Berne potential is short-ranged, most of the intermolecular 
interactions involve same-cell molecules; thus, only the relatively small number of molecules near cell boundaries will 
require interprocessor communications. The cell scheme and the required communications are somewhat difficult to 
implement, but provide very significant computational speedups. A computational scheme similar to ours is described 
in more detail in l40| (although we used a slightly different communication scheme in which an additional map tracking 
the specific subcells to be transferred between specific neighbors was implemented). We obtained timings virtually 
identical to those reported in the latter reference (on the order of 1 second per timestep with a 64 node partition of 
the Cray) . Our computation time increases linearly with the number of particles and with the number of processing 
elements, indicating good scalability of our code. 
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A preliminary step for locating defects is to break the system into a lattice of cubic bins. Note that the creation of 
a lattice is strictly for convenience in defect finding; the time evolution of the system allows for complete translational 
freedom. Even experiments, of course, have "binning" inherent in the resolution of the optical microscopy. Within 
each bin, the order parameter tensor Qap, Eq. (g), was calculated, its largest eigenvalue taken as the local order 
parameter S and the corresponding eigenvalue as the local director fi. The bin size was chosen so that the core size 
of the disclinations, determined from the distance over which S dropped significantly below the background value, 
was of the order of one lattice spacing. In our case, this resulted in a 16 x 16 x 8 lattice, each bin holding roughly 
30 molecules. It is with this lattice of orientation vectors that we began analyzing planes parallel to the x, y and z 
axes for the presence of defects. Note that while it is convenient to work with the orientation vectors on the lattice, 
we must remember that the actual directors are headless (expressing the symmetry upon rotation by 180° about an 
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axis perpendicular to the director) and so some care must be taken to account for this. 

A nice method of searching for discHnations was introduced in ||l^ (see also Q|). Consider the directors at the 
corners of a square (one of the faces of a cube) in our three-dimensional lattice. The idea is to track the course of 
the intersections of these vectors with the order parameter sphere (actually the projective plane RP2) as one moves 
around the corners of the real-space square (Fig. Eh. Starting with the intersection of Ha with the sphere, one then 
takes as the next point either the intersection of ns or — ns, whichever is closest to n^'s intersection. Once this point 
is determined, either he or —Ac is used, depending on the proximity to the previously defined point, and so on. Once 
the last point (from corner D) is determined, one looks at whether its intersection is in the same hemisphere as the 
starting point's. If so, no defect is present — the path in order parameter space is deformable to a single point, i.e., 
a uniform configuration. If the first and last points are in different hemispheres, however, then a disclination line is 
taken to cut through the center of the square and is oriented perpendicular to the plane of the square. 

To find escaped type-1 lines we performed a similar procedure, except that we measured the actual arclength 
swept out as one moves from each intersection to the next. Arclengths greater than n are counted as type-1 lines 
if the lattice square has not already been determined to hold a disclination (obviously, there is some overlap in the 
methods). The arclength tt corresponds to an escaped structure similar to that of Fig. g with an opening angle 
of about 30°. Experimentally, smaller opening angles are observable as type-1 lines, but smaller cutoff values of 
the arclength produced too many random type-1 line segments unconnected with each other or with disclinations, a 
situation not in accord with experiments. 




FIG. 4. The disclination-finding algorithm. The directors at the corners of the lattice cube face shown on the left are tracked 
on the order parameter space sphere shown on the right. The diameters AA',BB' ,CC' ,DD' correspond to the axes of the 
headless director at the real-space lattice sites A,B,C,D respectively. 



Finally, to look for monopoles, we used the method from ||^. Each of the six faces of a lattice cube is divided into 
two triangles and the corresponding directors are used to map out spherical triangles on the order parameter sphere. 
Total areas on the sphere greater than 2tt are considered to be monopoles and the defect is placed at the center of 
the corresponding lattice cube. Note that in all these defect-finding procedures, one must be careful to apply periodic 
boundary conditions to the edge lattice sites. 

One could also consider simulating the effect of crossed polarizers on individual planes. We used the method of 
|E3,E3[ which, phrased in the language of the Stokes parameters, uses Miiller matrices to simulate the effect of a 
group of molecules on the polarization of incoming light. In this method, one must set values for the ordinary and 
extraordinary refractive indices; we used typical experimental values between 1.5 and 2 
parameter, the ratio of the thickness of the cell to the wavelength of light, was chosen to be t 
calculated outgoing intensity for molecules oriented at 45° to the crossed polarizer directions equal to one; we used a 
value of 2.5. Visualizing the resulting contour plots is aided by choosing an exponential distribution of contour values 
in order to sharpen the dark areas (the "brushes" |4J]). This method did yield planes in which two clear brushes 
met at fairly well- localized points (Fig. 0), indicating the presence of a disclination but in general, the brushes and 
intersections were simply not well-determined enough to be useful. We estimate that an order of magnitude increase 
in the number of bins (corresponding to several million Gay-Berne particles) would be required to use this crossed 
polarizer approach quantitatively. 
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FIG. 5. Simulated crossed-polarizer image with actual director configuration superimposed. The crossed-polarizer image is 
the result of applying the Miiller matrix method to the single lattice plane shown. The disclination line (with topological charge 
1/2) at the top-center of the image is clearly indicated by two brushes. 
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A. Coarsening sequence 

With the methods described, we observed a coarsening sequence — compare Fig. |6| with similar figures in ||5[]"which 
exhibited most of the general behaviors observed experimentally B. An animation of our results is available on 
our web site |^^. Shortly after the quench, there was a dense tangle of defect lines. This tangle gradually thinned 
out and we could clearly identify and follow individual defect loops. With the exception of one wedge disclination 
line iQ] running through the sample, all of the disclination lines were of twist typejsee Fig. |j) and with periodic 
boundary conditions formed closed loops. The pre sence of twist lines is consistent [pSpJ] with the relative values of 
the elastic constants in the Gay-Berne nematic [Q, namely K2 < {Ki +K^)/2. Apart from the exception mentioned 
above we saw no evidence for dynamically generated wedge disclination lines that might contribute substantially to the 
structure factor. Combination, separation and collapse of the loops were all observed. The disclination loops appeared 
to experience minimal center-of-mass displacement and were relatively long-lived structures. Type-1 lines took the 
form of single line segments or small partial loops virtually always connected to disclination line segments and often 
forming bridges (much like the T intersections of J5|) between disclination segments from the same or distinct loops. 
Type-1 lines tended to fluctuate on much shorter time scales than disclinations, which seems reasonable given that 
the former are not topologically stable. One interesting observation is that type-1 lines often appeared as precursors 
to or remnants of the motion of disclination line segments. For example, the appearance of a type-1 line or several 
connected lines jutting out from a disclination was often followed by a kink or bend developing in the disclination 
line at that point. Similarly, the removal of kinks or bends often left behind type-1 lines for some period of time. 
The type-1 lines seemed to be initially defining and afterwards retaining a memory of the disclination path. Also, 
the emergence of distinct disclination loops from localized tangles often included the breaking of numerous type-1 
"bonds" between the loops. 
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FIG. 6. Coarsening sequence at times (a) i = 2, (b) t = 14, (c) t = 25, and (d) t = 37, with t = corresponding to the 
instantaneous temperature quench from the isotropic to the nematic phase. Filled squares represent point defects, thin lines 
represent type-1 lines and thick lines (for emphasis) represent disclinations. Note that with periodic boundary conditions, all 
disclination lines form closed loops. An animation of the coarsening sequence is available on our web site 1451. 
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FIG. 7. Indication of a twist disclination. A single lattice plane of directors is shown with disclinations indicated as thick 
lines. Dark areas indicate local directors perpendicular to the global director (along the vertical axis) while light areas indicate 
parallel orientations. The dark region in the center of the figure falls inside a disclination loop, clearly indicating a twist 
disclination. Compare with Fig. bib). 

The monopoles we observed fluctuated even more rapidly than the type-1 lines, although in many cases, their 
positions remained constant, on average, over relatively longer times. Because of their fluctuations, it is difficult 
to make any reliable statements about specific monopole behaviors such as monopole-antimonopole annihilation, for 
example. We never observed monopole formation upon disclination loop collapse, a result consistent with the presence 
of only twist disclination loops. All of the above processes are best observed in the animations we provide on our web 
site [43]. The total line length of disclination lines and type-1 lines as well as the number of monopoles is plotted in 
Fig. 8[as a function of time. The total number of disclination loops is plotted as a function of time in Fig. g. 
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FIG. 8. Time behavior of the various defect quantities: total length of disclination lines and type-1 lines and total number 
of monopoles. 
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FIG. 9. Number of disclination loops as a function of time. 



B. Real— space correlation function C{r,t) 



To calculate the correlation function C{r,t), Eq. (|3a|), we reduced our bin dimensions by a factor of 2 (yielding 
a 32 X 32 X 16 lattice size) to obtain a larger data set. We obtained curves for C{r, t) at times spanning the entire 
coarsening process. Motivated by the dynamical scaling hypothesis Eq. (|4^) we attempted to collapse our data to a 
single curve with appropriate rescaling of distances. From i = 13 (in units of thousands of steps since the temperature 
quench) until t — 40 when nearly all of the defects disappeared, the C{r, t) curves for different times collapse to 
a single curve (Fig. 11 (a)) upon rescaling distances by a length scale L{t) chosen so that C{r = L{t),t) = 1/2. 
This particular choice L{t) for the characteristic length scale was first suggested in |1^, and is the most accurate 
to implement numerically. The time dependence of L{t) is shown in Fig. |ll|(b). Our system is not large enough to 
extract a reliable power-law for the growth of L{t). However, according to the dynamic scaling hypothesis the length 
L{t) defined by the above criterion should differ at most by prefactors or subdominant contributions at late times 
from other characteristic length scales of the system. For example, as we noted in Sec. || the disclination line density 
should scale as (L(i))~^. In Fig. |l^ we plot the disclination line density as a function of L{t), over the range of times 
(t = 13 to i = 40) where we found good scaling of C{r, t). A least squares fit yields an exponent 1.99±0.23, consistent 
with dynamical scaling. However, the range of times over which coarsening occurs is too limited (due to the small 
system size) to allow us to fully assess the validity of the dynamical scaling hypothesis. Similarly, while the collapse 
of the correlation function data to a single scaling curve is consistent with the predictions of dynamical scaling, the 
range of distances and times is too limited to provide more definitive support for the hypothesis. Furthermore, as we 
discuss in the next section when we examine the structure factor dynamical scaling may in fact be breaking down in 
the range {r/L{t)) < 1, even though this is not evident from the scale of the plot of C{r, t). 
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FIG. 10. Log-log plot of disclination line density (line length per unit volume) versus the characteristic length L shown in 
Fig. [yib). The straight line is a least squares fit with slope —1.99 ±0.23. The range of L shown corresponds to the time range 
i = 13 to i = 40 where scaling behavior of C{r,t) is observed as shown in Fig. |ll|(a). 
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FIG. 11. a) Correlation function C{r,t) for times ranging from t = 12> to t = AQ with distances rescaled by the characteristic 
lengths defined by C{r — L(t)) = 1/2. (b) Time behavior of the characteristic length L{t). Note that the data in this figure 
was obtained by using a smaller bin size (corresponding to a 32 x 32 x 16 lattice) than the previous figures. 



C. Structure factor S{k,t) 



We computed the structure factor S{k,t), Eq.(3b), by first evaluating the Fourier transform of the nematic order 
parameter, Eq. (||): 



?a/3(k) = ^ E 2 



UiaUifi — -5ap 



exp(ik • r), 



(14) 



where V and N are the system volume and number of molecules respectively. As in Sec. VB we used a 32 x 32 x 16 
lattice. The wavevectors k have components which are multiples of the minimum values commensurate with the MD 
cell sizes in each direction. Motivated by the expression Eq. (4^) for the structure factor we plotted our results in 
log-log form. Representative results are shown in Figs. |2|, |l3 and |l^, where we have computed S{k) for values of 
k > 27r/16, the minimum commensurate value along the shortest dimension of the cell, and less than fc ~ 2. For 
values of k larger than 2 we are unable to fit our data to the long wavelength expression Eq. (3b) for S{k). Fig. n3 
corresponds to time t — 27 when there are still a sizable number of defects. Fig. |l^ corresponds to t = 40 near the 
end of the coarsening sequence, while Fig. U4 corresponds to i = 70 well beyond the end of the coarsening sequence. 
The data in the latter figure can be fit over nearly the entire range of fc by a power law iS'(fc) ^ fc~^-^, consistent 
with a purely thermal fluctuation contribution to the structure factor. On the other hand, we note that in the figure 
corresponding to a time midway through the coarsening process (Fig. (|l2|)), there is an apparent crossover in the 
behavior of S{k) as a function of k. For small values of fc, S'(fc) can be fit to the power law form, S'(fc) ~ fc"**-^, while 
at larger k the exponent is approximately 1.9. In Fig. nS we show the power-laws obtained at small and large values 
of fc during most of the coarsening sequence and beyond. This crossover behavior is consistent at least in part with 
the predictions of [Q. At large values of fc during the coarsening process (and at all values of fc when the process 
is complete), thermal fiuctuations dominate and we fit our data with an exponent close to 2. At small fc during the 
coarsening process our fit yields an exponent of approximately 4.5, clearly distinct from the thermal fluctuation form. 
As indicated in Fig. ||, the total disclination line length is in general an order of magnitude greater than the number of 
monopoles, i.e., pdisc — ^Opmonop- In spite of this order of magnitude difference in densities, we see from Eq. (pi) that 
the monopoles and disclinations should have comparable contributions to the structure factor in the range of fc values 
used in our plots. Thus, it is not clear why we obtain an exponent between 4 and 5, though one possibility is that 
we are seeing two-dimensional effects due to the anisotropic shape of our MD cell (recall that point disclinations in a 
two-dimensional nematic should yield a power law of 4). It is also possible that we have overestimated the number of 
monopoles, whose contribution to the structure factor would yield an exponent greater than 5. As discussed above in 
Sec. VA, the monopoles fluctuated quite rapidly, and it is possible that some apparent monopoles that we identified 



using the algorithm described in Sec. IV are not in fact topological defects 
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FIG. 12. Structure factor as a function of k at time t — 27. The exponents of the power-law fits at small and large k are 4.3 
and 1.8 respectively. At small k the structure factor is dominated by defects, while at large k thermal fluctuations dominate. 
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FIG. 13. Structure factor as a function of k at time t = 40, near the end of the coarsening sequence. The crossover between 
the thermal fluctuation regime at large k and the defect dominated regime at small k occurs at smaller values of k than at 
earlier times in the coarsening sequence (compare with Fig. |l2| ). Due to the relatively small number of data points in the defect 
dominated regime, we have not attempted a power-law fit. 
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FIG. 14. Structure factor as a function of k at time t = 70, after all of the defects have disappeared. The data is fit with an 
exponent of 1.7. With the coarsening process completed, only thermal fluctuations contribute to the structure factor. 
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FIG. 15. Values of the exponents used to fit the structure factor S{k) at small and large values of fc as a function of 
time. Topological defects dominate the small k behavior until the defects disappear around i = 40, while thermal fluctuations 
dominate at large k. At later times thermal fluctuations are the only contribution to S{k) for all k, and a single exponent flts 
the entire range of k. During the last stages of the coarsening sequence (between t — 30 and t — 40) we do not have sufficient 
numbers of data points to fit the small k behavior because the crossover to the thermal fiuctuation regime occurs at small 
k. However, the large k behavior continues to be fit well with an exponent of approximately 2 (see Fig. hS). The solid and 
dashed lines indicate the average values of the exponents used to fit the large and small k regimes of S{k) during the coarsening 
sequence; the values are 1.9 and 4.5 in these regimes respectively. The dashed-dotted curve indicates the average exponent, 
1.9, that fits all of the k data after the defects have disappeared. 



The crossover value of k separating the defect dominated and thermal fluctuation dominated regimes is of the order 
of magnitude predicted by Eq. (||), namely, k ^ {2'K^Kndisc/kBTY/^ ~ {'^^^TT^Pmono-plkBTY^^. This crossover value 
decreases with time, as can be seen by comparing Fig. |2| vifith the later time data of Fig. |l3[ The crossover value of k 
in the latter figure is about half of the corresponding value in the former figure, consistent with the relative densities 
of defects at the two times. 

As discussed in ref. [BOl we would expect the crossover to the thermal fiuctuation dominated regime to be accom- 
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panied by a breakdown of the dynamical scaling hypothesis because it assumes that thermal fluctuations play no role 
in the behavior of real-space or Fourier-space correlation functions. To test this expectation we used our data for 
S{k,t) during the time range spanning the coarsening process to plot the scaling function g defined in Eq. (4b). This 
plot is shown in Fig. |l6|, where we clearly see the breakdown of scaling for kL{t) greater than approximately 4 or 
5, corresponding to values of r/L(t) less than approximately one. Note that the numerical range of g is much larger 
than the range of /, the corresponding scaling function for C{r,t) (Eq. ([4a| ) and Fig. |ll|), so that the breakdown in 
scaling is easier to see in the structure factor data. The wider horizontal range for kL{t) compared to r/L{t) also 
makes the breakdown clearer. 
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FIG. 16. Plot of the dynamical scaling function <; = ^(fc, t)/(L(t))'' as a function of the scaling variable fcL(i), (see Eq. ([lb|) ) . 
Note the clear breakdown of scaling for large values of kL{t). 

Note in our plots, k > R~^, for nearly all values of the disclination loop radius R. Thus, in this regime we expect to 
see a power-law contribution to S{k) from the twist dislination loops. With a simulation of a larger system it might 
be possible to study the regime k <C R^^ where the twist disclination loops are expected to make no contribution to 
S{k), as discussed in Sec. ||. 
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In conclusion, we have shown that the Gay-Berne potential is fruitful for studying the behavior of the wide variety 
of topological defects generated in a quench from the isotropic to the nematic phase. At least for the Gay-Berne 
parameters chosen here, twist disclination loops were the dominant defects, and we did not, aside from possibly 
one isolated line, observe dynamically generated wedge disclination loops. This result, if it is not an artifact of our 
relatively small system size, has important implications for the interpretation of scattering experiments on quenched 
nematics, following upon the ideas of Ref. [Q. As we discussed in Sec. II measurements of the structure factor 
show scaling with an exponent between 5 and 6. Twist disclination lines yield an exponent of 5, whereas monopoles 
(which are relatively few in number in experimental systems) and wedge disclination lines (which are not energetically 
favorable) yield 6. Thus, it remains an open question as to why the exponent observed is greater than 5 over some 
appreciable range of wavevector. 

Our computed real-space correlation function exhibited good dynamical scaling over the limited range of distances 
available, though the structure factor appears to be provide a more sensitive test of scaling. In our structure factor 
data we could clearly see the breakdown of dynamical scaling and the crossover to the thermal fluctuation dominated 
behavior, in accord with the predictions of Zapotocky and Goldbart [pOJ . Clearly, simulations of even larger Gay-Berne 
systems would be of interest to further address the issues raised here. 
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